{
"cells": [
{
"cell_type": "markdown",
"id": "3d0a47e3-7c72-43fc-8563-d44cf886a24b",
"metadata": {},
"source": [
"# Introduction {#sec-zero}"
]
},
{
"cell_type": "markdown",
"id": "11700bf4",
"metadata": {},
"source": [
"---\n",
"format:\n",
" html:\n",
" other-links:\n",
" - text: This notebook\n",
" href: L1 Introduction to Math 5485.ipynb\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "24dc8e85",
"metadata": {},
"source": [
"\n",
"1. Syllabus + Canvas\n",
"2. ChimeIn\n",
"3. Jupyter notebooks\n",
"4. Basics of Numerical Analysis"
]
},
{
"cell_type": "markdown",
"id": "98d5dc48-aa86-4adf-81f4-6846acb068f3",
"metadata": {},
"source": [
"Jupyter Notebooks allow you to include Markdown, LaTex, and Julia (as well as python and other languages)"
]
},
{
"cell_type": "markdown",
"id": "fd40e03d-3a48-4f35-9c18-4fbbcc7c9bbd",
"metadata": {},
"source": [
"## Markdown"
]
},
{
"cell_type": "markdown",
"id": "e9092a73-504b-479d-aec5-528a2a6ab837",
"metadata": {},
"source": [
">Note: You may double click to see the syntax and edit cells. Press ctrl+enter to exit edit mode. \n",
">Pressing enter on the selected cell has the same effect as double clicking"
]
},
{
"cell_type": "markdown",
"id": "3b806ca9-3b86-4edc-8d52-c10a9fbd1141",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"Text in this cell can be **bold**, *italic*, ***bold italic***, quotes \n",
"\n",
">\"The only source of knowledge is experience\" \n",
">-- Einstein\n",
"\n",
"```\n",
"Preformatted text,\n",
" useful for\n",
"pseudocode\n",
"```\n",
"\n",
"Markdown cheat sheet: [link](https://www.markdownguide.org/basic-syntax/) "
]
},
{
"cell_type": "markdown",
"id": "6b065807-0ef0-4f0a-baa8-772e6212172c",
"metadata": {},
"source": [
"## LaTeX\n",
"\n",
"LaTeX allows you to display math: e.g. $\\pi \\approx \\frac{22}{7}$ and $$f(x) = f(a) + f'(a) (x - a) + \\frac{1}{2} f''(a) (x-a)^2 + R(x).$$\n",
"\n",
"LaTex cheat sheet: [link](https://quickref.me/latex.html) \n",
"Draw symbol to get latex syntax: [link](https://detexify.kirelabs.org/classify.html#google_vignette)"
]
},
{
"cell_type": "markdown",
"id": "dd55ba86-3acc-493f-8b1e-6a974e04dce6",
"metadata": {},
"source": [
"## Julia"
]
},
{
"cell_type": "markdown",
"id": "b8809e2e-19e8-49c8-9deb-b95a0d1fd3a2",
"metadata": {},
"source": [
">You can use greek symbols in Julia code: type \\pi and tab to get the symbol \n",
">You can show the result of a calculation by not ending a line with a semi-colon"
]
},
{
"cell_type": "code",
"execution_count": 60,
"id": "65def6a2-4c6e-47a2-bc77-56eb255296ca",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"π = 3.1415926535897..."
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"π"
]
},
{
"cell_type": "markdown",
"id": "b529b7c2-25b9-4b7b-80dc-2346f5b59667",
"metadata": {},
"source": [
">Note: you can add comments with #"
]
},
{
"cell_type": "code",
"execution_count": 61,
"id": "db3b6dab-12e1-4d45-9d6e-e96b8add1156",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.142857142857143"
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Define p\n",
"p = 22/7"
]
},
{
"cell_type": "markdown",
"id": "ad39b39b-358a-441d-a727-0f253ba7edd0",
"metadata": {},
"source": [
"Arrays start at 1:"
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "de3c03a6-ecf1-46fd-83ba-ad8558b29923",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×5 Matrix{Int64}:\n",
" 1 2 3 4 5"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = [1 2 3 4 5]"
]
},
{
"cell_type": "code",
"execution_count": 63,
"id": "c1b11be4-3601-4081-b163-4d1ee4c9798d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1"
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A[1]"
]
},
{
"cell_type": "markdown",
"id": "e8b2a769-303b-4eba-a708-0b99b1a53fdc",
"metadata": {},
"source": [
"Vectors:"
]
},
{
"cell_type": "code",
"execution_count": 64,
"id": "8b15d83c-aa52-46d2-89f4-a6b94130949c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Vector{Int64}:\n",
" 2\n",
" 4\n",
" 6\n",
" 8"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = [2, 4, 6, 8]"
]
},
{
"cell_type": "markdown",
"id": "8864e544-52c1-4392-a9ba-2c4351d39fdf",
"metadata": {},
"source": [
"Vector of objects:"
]
},
{
"cell_type": "code",
"execution_count": 65,
"id": "39b54b59-3d8e-44e8-aecb-cd80c6086dd6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Any}:\n",
" [1 2 … 4 5]\n",
" [2, 4, 6, 8]\n",
" \"hello world\""
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C = [ A , B , \"hello world\" ]"
]
},
{
"cell_type": "code",
"execution_count": 66,
"id": "17c2f160-509c-4bff-95b3-d912f9342abf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\"hello world\""
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"C[3]"
]
},
{
"cell_type": "markdown",
"id": "e7b52382-90e5-4a64-87c6-11dce4961656",
"metadata": {},
"source": [
"Matrices\n",
"\n",
"You can either specify the columns:"
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "59358e5f-3670-4369-a6a5-d392434c6b7f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Int64}:\n",
" 1 2\n",
" 4 5"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = [ [1, 4] [2, 5] ]"
]
},
{
"cell_type": "markdown",
"id": "e396fedc-a220-442a-9eff-51f627e2dcc1",
"metadata": {},
"source": [
"Or the rows: "
]
},
{
"cell_type": "code",
"execution_count": 68,
"id": "9829675a-2f68-4188-a5fe-86b9f9e66dbf",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Int64}:\n",
" 1 2\n",
" 4 5"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Y = [ 1 2 ; 4 5 ]"
]
},
{
"cell_type": "markdown",
"id": "d32edf11-d4a2-4eef-b86a-6feeb1b2aa84",
"metadata": {},
"source": [
"Linear Algebra:"
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "e248b5ef-b751-4950-9895-92bde3f9214d",
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra"
]
},
{
"cell_type": "code",
"execution_count": 70,
"id": "01aec85f-a449-425f-bdb3-a4bb483414a1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-3.0"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"det( X )"
]
},
{
"cell_type": "code",
"execution_count": 71,
"id": "0b34234c-de70-4f93-bc78-b576daef4e0a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tr( X )"
]
},
{
"cell_type": "code",
"execution_count": 72,
"id": "3a5153e3-f58e-45f0-b13f-ef246cc369ca",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Float64}:\n",
" -1.66667 0.666667\n",
" 1.33333 -0.333333"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(-1/3) * [ 5 -2 ; -4 1 ]"
]
},
{
"cell_type": "code",
"execution_count": 73,
"id": "9795e903-9c8b-4b21-aba4-954adaef9329",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Float64}:\n",
" -1.66667 0.666667\n",
" 1.33333 -0.333333"
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"inv( X )"
]
},
{
"cell_type": "code",
"execution_count": 74,
"id": "52357743-8092-480d-86b0-fa49f2c18393",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2-element Vector{Float64}:\n",
" -0.4641016151377544\n",
" 6.464101615137754"
]
},
"execution_count": 74,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ϵ = eigvals( X )"
]
},
{
"cell_type": "code",
"execution_count": 75,
"id": "d879a0f5-9298-40b2-8c61-bb126eba299f",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×2 Matrix{Float64}:\n",
" -0.806898 -0.343724\n",
" 0.59069 -0.939071"
]
},
"execution_count": 75,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v = eigvecs( X )"
]
},
{
"cell_type": "code",
"execution_count": 76,
"id": "1ac4a1ea-1273-4cc9-bb88-605217e62980",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.879951945339055e-16"
]
},
"execution_count": 76,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"v1 = v[1:2,1]; # gets the first column of v\n",
"norm( ( X - ϵ[1] * Diagonal([1, 1]) ) * v1 )"
]
},
{
"cell_type": "markdown",
"id": "f4082477-ecfd-4e4c-aee4-77e8cad66978",
"metadata": {},
"source": [
"Functions"
]
},
{
"cell_type": "code",
"execution_count": 77,
"id": "2003c7a0-68a7-46ec-b6e9-7a12a9cc8db8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-1.0"
]
},
"execution_count": 77,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f = x -> cos( x );\n",
"f( π )"
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "b3519246-2692-4513-873d-a3dc8192882b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" -1.0\n",
" 1.0\n",
" -1.0"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# vectorisation = entry-wise functuions --> use dot!\n",
"f.( [-π, 0, π] )"
]
},
{
"cell_type": "code",
"execution_count": 79,
"id": "6b56b002-fb59-41fa-9fc4-55f1e54cf0dd",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3-element Vector{Float64}:\n",
" -1.0\n",
" 1.0\n",
" -1.0"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# you can also specify the vectorisation at the start:\n",
"@. f( [-π, 0, π] )"
]
},
{
"cell_type": "markdown",
"id": "24d6a0ba-c1f2-402e-8caa-137abccc175f",
"metadata": {},
"source": [
"Plotting functions:"
]
},
{
"cell_type": "code",
"execution_count": 80,
"id": "a03ae8ef-8ccb-455a-9e3c-34704dfc422c",
"metadata": {},
"outputs": [],
"source": [
"# We need Plots to plot graphs\n",
"using Plots\n",
"# LaTeXStrings lets you put latex in graph titles and labels\n",
"using LaTeXStrings"
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "a19628b1-df4e-4ee9-9ea4-1f64c7c3c00d",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = -5:.01:1;\n",
"Y = X .* exp.(X) ;\n",
" \n",
"plot( X, Y, label=L\"f\", xlabel=L\"w\", title=L\"f(w) =w e^w\") \n",
"\n",
"# L\"\" is a LaTeX string\n",
"# \"\" is a normal srtring"
]
},
{
"cell_type": "markdown",
"id": "452b1522-7ec9-47cc-a1cc-e49f8bd12b90",
"metadata": {},
"source": [
"Element-wise operations can also be written using the following:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bf0ece6c-77c7-4702-9090-3ddf8d8b5495",
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#| label: fig-1\n",
"#| fig-cap: \"Ploting functions\"\n",
"#| column: margin\n",
"\n",
"X = range(0, 10, length=100); \n",
"Y1 = @. sin(X) / X; # this is the same as Y1 = sin.(X) ./ X; \n",
"Y2 = cos.(X); \n",
"plot( X, [Y1 Y2], xlabel=L\"x\", label=[L\"\\frac{\\sin x}{x}\" L\"\\cos x\"], linewidth=[4 2], lc=[:green :blue]) \n",
"\n",
"\n",
"# this the same as\n",
"plot( X, Y1, xlabel=L\"x\", label=L\"\\frac{\\sin x}{x}\", linewidth=4, lc=:green) \n",
"plot!( X, Y2, label=L\"\\cos x\", linewidth=2, lc=:blue) \n",
"# without the !, only the most recent function will be plotted\n"
]
},
{
"cell_type": "code",
"execution_count": 83,
"id": "fec8b9f9-b079-44e3-bd10-853b83f78d56",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f1(x) = 1 / (1 + x^2);\n",
"f2(x) = 1 / (1 + sin(x)^2);\n",
"plot(f1, -π, π, lw=3, label = L\"f_1(x) = \\frac{1}{1+x^2}\");\n",
"plot!(f2, -π, π, lw = 3, label = L\"f_2(x) = \\frac{1}{1+\\sin^2(x)}\", \n",
" size = (600, 250), legend = :outertopright )"
]
},
{
"cell_type": "markdown",
"id": "4579e5b3-3dbf-459b-997e-6a36e819d04c",
"metadata": {},
"source": [
"# Chapter 1: Basics of Numerical Analysis {#sec-one}\n",
"\n",
"- Absolute and relative error\n",
"- Rate of convergence\n",
"- Using plots to show rate of convergence"
]
},
{
"cell_type": "markdown",
"id": "291ba2d9-59df-48cc-9fb1-35c1f3fc0cd7",
"metadata": {},
"source": [
"Suppose we are approximating $x$ with $\\widetilde{x}$. Here, $x$ could be e.g. a root (or fixed point) of a function, the value of an integral, the solution of an ODE or PDE evaluated at a point, .... . "
]
},
{
"cell_type": "markdown",
"id": "edf22d61",
"metadata": {},
"source": [
"> **Definition 1 (Errors).** We define $$\\text{(Actual) Error:} \\qquad x - \\widetilde{x},$$ $$\\text{Absolute Error:} \\qquad |x - \\widetilde{x}|,$$ and, when $x \\not= 0$, we define $$\\text{Relative Error:} \\qquad \\frac{|x - \\widetilde{x}|}{|x|}$$"
]
},
{
"cell_type": "markdown",
"id": "22b4689f",
"metadata": {},
"source": [
"> **Example 1.** *(i)* If $x = 1.5$ and $\\widetilde{x} = 1$, then the absolute error is $0.5$ and the relative error is $\\frac{0.5}{1.5} = \\frac13$. \n",
">\n",
"> *(ii)* If $x = 1.5 \\times 10^{42}$ and $\\widetilde{x} = 10^{42}$, then the absolute error is $0.5 \\times 10^{42}$ and the relative error is $\\frac13$.\n",
">\n",
"> *(iii)* If $x = 1.5 \\times 10^{-42}$ and $\\widetilde{x} = 10^{-42}$, then the absolute error is $0.5 \\times 10^{-42}$ and the relative error is $\\frac13$."
]
},
{
"cell_type": "markdown",
"id": "4ebe53fe",
"metadata": {},
"source": [
"The relative error is a better measure for small or large $x$. e.g. approximating Planck's constant $6.62607015 \\times 10^{-34} \\text{m}^2 \\text{kg } \\text{s}^{-1}$ with $0 \\text{m}^2 \\text{kg } \\text{s}^{-1}$ would give an absolute error of $\\approx 6 \\times 10^{-34}$ (which seems small) but a relative error of $1$. In practice $0$ is a bad approximation to Planck's constant!"
]
},
{
"cell_type": "markdown",
"id": "ddb6e8b2",
"metadata": {},
"source": [
">**Definition 2 (Number of accurate digits)**. Suppose $x \\not= 0$ is approximated by $\\widetilde{x}$. Then, we say this approximation is accurate to \n",
">$$\\eta := -\\log_{10} \\left| \\frac{x-\\widetilde{x}}{x} \\right| $$\n",
">digits."
]
},
{
"cell_type": "markdown",
"id": "fbe8492c",
"metadata": {},
"source": [
">**Example 2.** Approximating $e := \\lim_{n\\to\\infty} \\big( 1 + \\frac1n \\big)^n \\approx 2.7182818.....$ with $3, 2.7,$ and $2.718$ gives $\\eta = 0.98, 2.17$, and $3.98$ (to $2$ decimal places), respectively."
]
},
{
"cell_type": "markdown",
"id": "884d8567",
"metadata": {},
"source": [
"We now consider a sequence of approximations $x_n$ converging to $x$. That is, for all $\\varepsilon > 0$, there exists $N=N_\\varepsilon$ such that $|x - x_n| < \\varepsilon$ for all $n > N$. In numerical analysis, we are interested in *(i)* whether our numerical schemes converge, and *(ii)* the rate of convergence: "
]
},
{
"cell_type": "markdown",
"id": "c6a996b0",
"metadata": {},
"source": [
">**Definition 2 (Rate of convergence).** Suppose $\\varepsilon_n \\geq 0$ and $\\varepsilon_n \\to 0$. Then, $x_n \\to x$ with rate $O(\\varepsilon_n)$ if there exists $C, N$ such that $|x - x_n| \\leq C \\varepsilon_n$ for all $n > N$. We also write $x_n = x + O(\\varepsilon_n)$."
]
},
{
"cell_type": "markdown",
"id": "1648de06",
"metadata": {},
"source": [
">**Exercise 1.** Show that $\\frac{n + \\cos(n)}{n^2} = O(\\frac1n)$ and $\\frac{n^2+3}{n^{42}} = O(\\frac1{n^{40}})$\n"
]
},
{
"cell_type": "markdown",
"id": "541e06ae",
"metadata": {},
"source": [
"When plotting graphs, $O(n^{-k})$ decay (that is, **algebraic** decay) can be seen as straight lines on log-log plots ($-k$ is the slope in such graphs):"
]
},
{
"cell_type": "code",
"execution_count": 84,
"id": "7f0767c3",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = 1:100;\n",
"plot( N.^(-1) , xlabel=\"n\", xaxis=:log, yaxis=:log, label = L\"n^{-1}\")\n",
"plot!( N.^(-2), label = L\"n^{-2}\")\n",
"plot!( N.^(-5), label = L\"n^{-5}\")"
]
},
{
"cell_type": "markdown",
"id": "89fdfc1b",
"metadata": {},
"source": [
"On the other hand, $O(e^{-c N})$ (that is, **exponential** decay) can be seen as straight lines on a semi-log y plot (and $-c$ is the gradient):"
]
},
{
"cell_type": "code",
"execution_count": 85,
"id": "d9dcce18",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot( exp.(-N) , xlabel=\"n\", yaxis=:log, label = L\"e^{-n}\")\n",
"plot!( exp.(-2 .* N ), label = L\"e^{-2n}\")\n",
"plot!( exp.(-5 .* N ) , label = L\"e^{-5n}\")"
]
},
{
"cell_type": "markdown",
"id": "be06175b",
"metadata": {},
"source": [
"# Exercises"
]
},
{
"cell_type": "markdown",
"id": "1bb51f88-2429-4849-b133-043e54d6ef6b",
"metadata": {},
"source": [
"1. $e_N := \\big( 1 + \\frac{1}{N} \\big)^N$ is an approximation to Euler's number $e := \\exp(1)$. What is the error, absolute error, and relative error in approximating $e$ with $e_N$ for $N = 1,2,...,10$? What are the number of accurate digits in this approximation?\n",
"2. What is the rate of convergence of this method?"
]
},
{
"cell_type": "code",
"execution_count": 86,
"id": "89b5c81f-7555-48af-9164-68c7062957ba",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2.718281828459045"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e = exp(1)"
]
},
{
"cell_type": "code",
"execution_count": 87,
"id": "57f26033-3822-4550-ab63-413bb81b91df",
"metadata": {},
"outputs": [],
"source": [
"N = 1:100;\n",
"err = @. abs( exp(1) - (1 + 1/N)^N );"
]
},
{
"cell_type": "markdown",
"id": "1015d55f-acde-4a9a-a38b-14e549dca05a",
"metadata": {},
"source": [
"Error = absolute error (in this case)"
]
},
{
"cell_type": "code",
"execution_count": 88,
"id": "14c4be48-e88f-495f-b033-53db5d65e4fb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10-element Vector{Float64}:\n",
" 0.7182818284590451\n",
" 0.4682818284590451\n",
" 0.34791145808867485\n",
" 0.2768755784590451\n",
" 0.22996182845904567\n",
" 0.19665545671693163\n",
" 0.17178213141833298\n",
" 0.1524973145086972\n",
" 0.13710703674584668\n",
" 0.12453936835904278"
]
},
"execution_count": 88,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"err[1:10]"
]
},
{
"cell_type": "markdown",
"id": "4f67f8f2-7d3e-4c58-a0b3-9b64f71c0eac",
"metadata": {},
"source": [
"Relative error"
]
},
{
"cell_type": "code",
"execution_count": 89,
"id": "92c0058c-951e-4bd3-9f29-02181c43ad64",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10-element Vector{Float64}:\n",
" 0.26424111765711533\n",
" 0.17227125736425472\n",
" 0.12798947277880338\n",
" 0.10185683307753335\n",
" 0.08459822894427681\n",
" 0.07234549952033957\n",
" 0.0631951145094156\n",
" 0.05610062684160521\n",
" 0.050438860058734485\n",
" 0.045815473235769066"
]
},
"execution_count": 89,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"err[1:10]/exp(1)"
]
},
{
"cell_type": "markdown",
"id": "d98694e0-7bbd-464b-8a75-fa9bb922e615",
"metadata": {},
"source": [
"Number of correct digits:"
]
},
{
"cell_type": "code",
"execution_count": 90,
"id": "ea9d6daf-c9d7-4a7c-a611-49c8e8f36128",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"10-element Vector{Float64}:\n",
" 0.5779996023832055\n",
" 0.7637871764659977\n",
" 0.8928257498997241\n",
" 0.992009830990522\n",
" 1.072638728778693\n",
" 1.1405884803607715\n",
" 1.199316494876063\n",
" 1.2510322861176486\n",
" 1.297234737241614\n",
" 1.338987823174671"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"η = - log10.( err[1:10]/exp(1) )"
]
},
{
"cell_type": "markdown",
"id": "93085666-cd61-41f9-9a57-66cc778a39b5",
"metadata": {},
"source": [
"Rate of convergence:"
]
},
{
"cell_type": "code",
"execution_count": 91,
"id": "31dc7b83-f73b-4cdf-b866-a01433e6dbf8",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot( err , xlabel=\"N\", title=\"absolute error\", legend=false)\n",
"plot!( N.^(-1) )"
]
},
{
"cell_type": "code",
"execution_count": 92,
"id": "a09b8d12-5b05-4544-b5fb-e6343e71207f",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot( err , xaxis=:log, yaxis=:log, xlabel=\"N\", title=\"absolute error\", legend=false)\n",
"plot!( N.^(-1) )"
]
},
{
"cell_type": "markdown",
"id": "a380a17a-671f-4bba-9484-cc2cbed21344",
"metadata": {},
"source": [
"3. Let $\\gamma_n := \\sum_{k=1}^n \\frac{1}{k} - \\log n$ and define the Euler–Mascheroni constant as $\\gamma := \\lim_{n\\to\\infty} \\gamma_n$. What is the absolute error in approximating $\\gamma$ with $\\gamma_{1000}$?\n",
"4. What is the rate of convergence of this approximation to $\\gamma$?"
]
},
{
"cell_type": "code",
"execution_count": 93,
"id": "8174bb48-9065-4d83-b40d-859476fa6d8c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"γ = 0.5772156649015..."
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"γ = Base.MathConstants.eulergamma"
]
},
{
"cell_type": "code",
"execution_count": 94,
"id": "4d5d5c46-d79f-43c1-b147-2aa22a1a5db5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"gamma[1000] = 0.5777155815682078\n",
"Absolute error = 0.0004999166666749266\n"
]
}
],
"source": [
"Nmax = 1000;\n",
"\n",
"N = 1:Nmax;\n",
"γN = zeros( Nmax );\n",
"γN[1] = 1;\n",
"for n in 2:Nmax\n",
" γN[n] = (1/n) + γN[n-1] + log( (n-1)/n );\n",
"end\n",
"\n",
"println( \"gamma[1000] = \", γN[1000] )\n",
"println( \"Absolute error = \" , abs( γ - γN[1000] ) )"
]
},
{
"cell_type": "markdown",
"id": "9ae9d733-05c3-434c-90e1-54915076d7d9",
"metadata": {},
"source": [
"Number of accurate digits:"
]
},
{
"cell_type": "code",
"execution_count": 95,
"id": "0dd2f42c-de2d-4fb8-8958-e061d121ef97",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3.0624404928861617"
]
},
"execution_count": 95,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"η = - log10.( abs( γ - γN[1000] )/ γ ) "
]
},
{
"cell_type": "code",
"execution_count": 96,
"id": "5de7f1bb-0a32-476e-977c-36cc2f54693e",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot( γN .- γ , xaxis=:log, yaxis=:log , legend = false, lw = 2) \n",
"plot!( N.^(-1) , linestyle=:dash)"
]
},
{
"cell_type": "markdown",
"id": "ea8fd20a",
"metadata": {},
"source": [
"5. [This was just so that we can remember some calculus and also rigorously prove the rate of convergence that we have seen above using a rigorous error estimate] \n",
" Prove that $|\\gamma - \\gamma_n| \\leq \\frac1n$."
]
},
{
"cell_type": "markdown",
"id": "dfae1e3a",
"metadata": {},
"source": [
"*Proof.* First, notice that $\n",
"\\begin{align} \n",
" \\gamma_n - \\gamma_N \n",
" %\n",
" &= \\log \\frac Nn - \\sum_{k=n+1}^N \\frac{1}{k}\\\\\n",
" %\n",
" &= \\sum_{k=n+1}^N \\left( \\log \\frac{k}{k-1} - \\frac{1}{k} \\right).\n",
"\\end{align}$\n",
"Here, $\\frac1k$ is the area under the rectangle with base $[k-1,k]$ and height $\\frac1k$ and $\\log \\frac{k}{k-1}$ is the integral of $\\frac1x$ between $x = k-1$ and $x = k$: "
]
},
{
"cell_type": "code",
"execution_count": 97,
"id": "0794b720",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 97,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k=3;\n",
"x = range(k-1,k, 100);\n",
"plot( x, [x.^(-1) ones(100)./k (x .+ 1).^(-1)], \n",
" ylims=(0,1/(k-1)), \n",
" label=[L\"y=\\frac{1}{x}\" L\"\\frac{1}{k}\" L\"y=\\frac{1}{x+1}\"] )"
]
},
{
"cell_type": "markdown",
"id": "0d42bcf3",
"metadata": {},
"source": [
"In the graph above, $\\frac1k$ is the area under the red line, $\\log \\frac{k}{k-1}$ is the area under the blue line. We upper bound the absolute value of the error between these areas, by computing the area between $y = \\frac1x$ and $y=\\frac1{x+1}$:\n",
"\n",
"\\begin{align}\n",
" 0 \\leq \\gamma_n - \\gamma_N \n",
" %\n",
" &\\leq \\sum_{k=n+1}^N \\left( \\log \\frac{k}{k-1} - \\frac{1}{k} \\right) \\\\\n",
" %\n",
" &\\leq \\sum_{k=n+1}^N \\int_{k-1}^k \\left( \\frac{1}{x} - \\frac{1}{x + 1} \\right) \\mathrm{d}x \n",
" %\n",
" = \\int_{n}^\\infty \\frac{1}{x(x+1)} \\mathrm{d}x \\\\\n",
" %\n",
" &\\leq \\int_{n}^\\infty \\frac{1}{x^2} \\mathrm{d}x = \\frac1n\n",
"\\end{align}\n",
"\n",
"Here, we have shown that $(\\gamma_n)$ is a Cauchy sequence and thus it converges, and we have the same estimate for $\\gamma_n - \\gamma$ by taking the limit in $N$. $\\text{ }\\square$ "
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.7.1",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 5
}